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Abstract 

A recently developed analytical method for systematic improvement of the con- 
vergence of path integrals is used to derive a generalization of Euler's summation 
formula for path integrals. The first p terms in this formula improve convergence of 
path integrals to the continuum limit from 1/N to l/N^, where N is the coarseness 
of the discretization. Monte Carlo simulations performed on several different models 
show that the analytically derived speedup holds. 
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1 Introduction 



Feynman's path integrals [1,2] represent a compact and rich formalism for 
dealing with quantum theories. They have proved to be powerful tools for 
investigating symmetries, deriving non-perturbative results, delineating con- 
nections between different theories and different sectors of theories. Their flex- 
ibility and intuitive appeal have allowed us to generalize quantization to ever 
more complicated systems and have led to a rich cross fertilization of ideas 
between high energy and condensed matter physics [3,4]. Today, they are used 
both analytically and numerically [5,6,7,8,9] in many other areas of physics, 
in chemistry and materials science, as well as in quantitative finance [10]. 
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Unfortunately we still have very little knowledge of the precise mathematical 
properties of path integrals. As a result only an extremely limited number 
of path integrals can be solved exactly. Although the functional formalism 
has been instrumental for deriving many general approximation techniques, 
along with a host of model-specific approximations, there remain many mod- 
els of interest that need to be treated numerically (e.g. using Monte Carlo 
simulations) . 

Numerical integrations of path integrals have, however, proven to be notori- 
ously demanding of computing time. For this reason several research groups 
have worked at improving the convergence of path integrals. Until recently 
the best available result for partition functions of a generic A^-fold discretized 
theory led to a convergence [11,12,13]. A new systematic analysis of the 
relation of discretizations of different coarseness [14,15] has made it possible 
to dramatically improve convergence for general transition amplitudes, not 
just partition functions. A result of this investigation has been a procedure 
for constructing a series of effective actions S^^^ having the same continuum 
limit as the starting action S, but which approach that limit as l/N^. Explicit 
expressions for these effective actions have so far been constructed for p < 9 
[15,16]. In the current paper we will build on this derivation (and simplify it to 
some extent) and will cast the new analytical input in the form of a generalized 
Euler summation formula for path integrals. It is our belief that the existence 
of such a general formula strongly hints at the possibility of (currently unseen) 
simplifications that might make it possible to set up a rigorous theory of path 
integration. In section 2 we present a simple derivation of Euler 's summation 
formula for ordinary integrals as a useful guide to the generalization to path 
integrals that is given in sections 3 and 4. 



2 Euler summation formula for ordiuciry integrals 

The current status of the development of the path integral formalism is quite 
similar to that of ordinary integrals before the setting up of integration theory 
by Riemmann. In those days integrals were calculated directly from the defin- 
ing formula, i.e. one looked at a specific discretization of the integral (Darboux 
sum), attempted to do the sum explicitly, and finally tried to calculate the 
continuum limit. For example, 

^ JV 

/[/] ^ / f{t)dt = hm IM , where /,v[/] = E f{tn)eN , (1) 



En = T/N and tn = ne^- It goes without saying that done this way, even the 
simplest ordinary integrals presented a challenge. The mathematicians of the 
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18th century did not have computers at their disposal or the development of in- 
tegration theory might have come much later, i.e. they might have succumbed 
to doing brute force numerical calculations of integrals of all but the simplest 
functions. The problem with these hypothetical numerical calculations would 
have been two fold: they would have been inefficient (the discretized sums 
converge slowly to the continuum value), and they would have worked (thus 
quite probably slowing down the further development of integration theory). 
Luckily, this early numerical road was not open. The last great step in the 
development of integration before Riemmann was made by Euler. 

Discretization is not unique. This makes it possible to change f{t) to some 
other function (adding terms proportional to ejv, e%, etc.) without changing 
the integral. Let us assume that f*(t) is such an equivalent function with the 
added property that the sums InII*] do not depend on N. In fact we shall 
present a way of exphcitly constructing f*{t) for any given f{t). We first look 
at the simple case oi f{t) — 1. Now 

lN[l]^J2eM^T , (2) 

n=l 



which is already A^-indepcndcnt. Hence, in this case, all the additional terms 
vanish. Note that f*{t) is completely determined by the original function f{t) 
(and by ejv), so that the additional terms necessarily depend only on the 
derivatives /', /", etc. 

The second step is to take f{t) = t. In this case we get 
^ ^ ^ N(N + 1) 



From this it follows that lN[t ~ ^] = "y- Therefore, up to /" and higher 
derivatives of / that all vanish for linear f{t), we have f*{t) — f{t) — ^f'{t). 

We continue this procedure by looking at f{t) —t^. In this case we find 
r r.2i >^.2 N{N + 1){2N +1) 



It follows that /Ar[t^ — epitn — = In terms of /* this gives f*{t) = 

fit) — ^f'(t) ^f"{t) + . . .. The additional terms now depend on higher 

powers of as well as on higher derivatives and are determined by considering 
/Ar[t^], and so on. In this way we have constructed a procedure for finding f*{t) 
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for any given f{t). Remembering that InU*] does not depend on N we find 
/ mdt = E fiQe^ - f ^ n^er, - ^ E f"iin)eN + ■ ■ ■ ■ (5) 

Q n=l ^ n=l n=l 

This is the well-known Euler summation formula. We may also write it more 
compactly as 

m-lN[f^^]+0{el), (6) 

where f^^ is the truncation of /* to the first p terms. The Euler formula gives 
the analytical relation between integrals and their discretized sums. Looked 
at numerically, this formula allows us to increase the speed of convergence of 
discretized expressions to the continuum limit. For example, in the defining 
relation the discretized expressions differ from the continuum by a term of 
order 0{1/N). By using the Euler sum formula with p terms we can reduce 
that error to 0(1/ N^). All that is needed to do this is that the integrand is 
difi^crcntiable p — 1 times, the following sections we will generalize the above 
approach to path integrals. 



3 General properties of path integrals 



In the functional formalism the quantum mechanical amplitude A{a, 6; T) = 
(6|e^^^|a) is given in terms of a path integral which is simply the N ^ oo 
limit of the {N — 1)- fold integral expression 



AN{a,b;T) = {^—] j dqi ■ ■ ■ dgN-ie'^'^ . (7) 



The Euclidean time interval [0, T] has been subdivided into equal time steps 
of length = T/N, with qq — a and qn — b. Sn is the naively discretized 
action of the theory. We focus on actions of the form 

^ 1 

^2 







dt[-q' + V{q)) , (8) 



whose naive discretization is simply 



N-l / ^2 



n=0 
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where 5„ = - g„, K = V{qn), and g„ = |(?n+i + ?n)- We use units in 
which h and particle mass equal 1. 

As was the case with ordinary integrals the definition of the path integrals 
also makes it necessary to make the transition from the continuum to the dis- 
cretized theory, a process that is far from unique. For theories described by 
eq. (8) we have the freedom to choose any point in in which to eval- 

uate the potential without changing physics - the discretized amplitudes do 
differ, but they tend to the same continuum hmit. The calculations we present 
turn out to be simplest in the mid-point prescription where the potential V 
is evaluated at A more important freedom related to our choice of dis- 
cretized action has to do with the possibility of introducing additional terms 
that explicitly vanish in the continuum limit. Actions with such additional 
terms will be called effective. For example, the term Y.n=o '^n SndiQn): where 
g is regular when eAr — > 0, does not change the continuum physics since it 
goes over into Jo dtq'^ g{q), i.e. it vanishes as e%. Such terms do not change 
the physics, but they do affect the speed of convergence. A systematic study 
of the relation between different discretizations of the same path integral will 
allow us to explicitly construct a series of effective actions with progressively 
faster convergence to the continuum. Before we do this we will parallel the 
derivation in the previous section and derive some general properties of the 
best effective action. 

The amplitude A{a,b;T) of some theory with action S satisfies 

A{a, b;T) = / dqi--- dqn-iA{b, q^-i, e^) ■ ■ ■ A{qi, a; ejv) , (10) 



for all A^. This general relation is a direct consequence of the linearity of states 
in a quantum theory. In analogy with ordinary integrals let us now suppose 
that there exists an effective action S* that is equivalent to S (i.e that leads to 
the same continuum limit for all path integrals) with the additional property 
that its A^-fold discretized amplitude A*j^{a,b]T) does not depend on N, i.e. 
that satisfies 

A*j,{a,b-T) = A{a,b-T) . (11) 



As was the case in the previous section we will in fact construct a general 
procedure for evaluating this effective action. For actions of the form given in 
eq. (8) we may write the amplitude as 



^(?n+i, qn] ejv) = (2^) (~2e^) "^^^"+^' ' ^^'^^ 
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where the reduced amphtude ^— >laseAr— >0. Writing S^f as 

N-l / .2 \ 

Sn=j: {^+^nw:] , (13) 

n=o V^^AT J 

and using eq. (7), (10) and (11) we find 

exp (-ejvW*) = A(qn+i, Qn] ejv) ■ (14) 

Note that W* is reminiscent of some effective potential, so it should depend 
on g„, however, from the above relation we see that it must also depend on 
6n- In addition, W* also has an explicit dependence on the discrete time step 
En, hence 

W:^W*{5n,qn;eN) . (15) 



As we have seen, the above functional form is a direct consequence of the 
linearity of quantum theory. The equivalence of S and S* implies that W* 
V{q) when and 5 go to zero. The final general property of W* follows from 
the reality of amplitudes in the Euclidean formalism. Using the hermiticity of 
the Hamiltonian we find A{a, b; T) = A{a, b; T)t = (6|e-^^|a)t = (a|e-^^|6) = 
A{b, a; T). In terms of W* this gives us 

W*{5n, Qn, ejv) = W*{-dn, Qn, Cjv) , (16) 

or, said another way, only even powers of Sn are present in the expansion of 
W*: 

W*{Sn, qn, €n) = goiqn, ^n) + ^iQli^n] ^n) + 5n92{qn] Cjv) + ■ ■ ■ ■ (17) 



All the functions are regular in the e — > limit. The link to the starting 
theory is now simply go{qn', ^n) — V{qn) as goes to zero. This concludes 
the general properties of W*. The remaining properties will be analyzed in 
the following section by studying the relation of discretizations of different 
coarseness. 



4 Euler summation formula for path integrals 



We start by studying the relation between the 2iV-fold and iV-fold discretiza- 
tions of the same theory. Prom eq. (7) we see that we can write the 2A^-fold 
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amplitude as an A'"-fold amplitude given in terms of a new action Sn deter- 
mined by 



e 



2 ^ 

-^^ = (^) ' J d^i--- e"^2~ , (18) 



where S2N is the 2A'"-fold discretization of the starting action. We have written 
the 2A^-fold discrctizcd coordinates Qo, Qi-i ■ ■ ■ , Q2N in terms of g's and x's in 
the following way: Q2k = Qk and Q2k~-i = ^k- Note that we have go = ct, <1n = b, 
while the N — 1 remaining g's play the role of the dynamical coordinates in the 
/^-fold discretized theory. The x's are the N remaining intermediate points 
that we integrate over in eq. (18). It is not difficult to see that if we use the 
naively discretized action Sn one obtains for Sn an expression that is not of 
the same form as S'^v- 

Having in mind the results of the previous section it is best to use the effective 
action 

N-l / p \ 

S*N=T.{7r^ + ^NW*iS^, er,) , (19) 



which gives the same result for both the 2iV-fold and iV-fold discretizations. 
Therefore, in this case we get 



2 ~ 

e-s*^ ^ ( ) ' I dxi---dxN e~^2^ . (20) 

Xttcn J J 



From this one easily finds 

1 +00 

— oo 



where 



- — lnF(x;eiv) = QqI^ ;ejvj +^o(^ — ^ — ;eivj + 

+ (q„+i - xf gi (^^!^±11j^: ; + {x - qnf Qi {^-^-^ ; + . . . (22) 



The above integral equation can be solved for the simple cases of a free particle 
and a harmonic oscillator, and gives the well known results. Note however 
that for a general case the integral in eq. (21) is in a form that is ideal for 
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an asymptotic expansion [17]. The time step is playing the role of small 
parameter (in complete parallel to the role h plays in standard semi-classical, 
or loop, expansion). As is usual, the above asymptotic expansion is carried 
through by first Taylor expanding F (qn + y, ^) around Qn and then by doing 
the remaining Gaussian integrals. Assuming that < 1 (i.e. N > T) we have 



= In 



E 

m=0 



F^'™) {qn, f 
(2m)!! 
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(23) 



Note that F^'^'^\x;eN) denotes the corresponding derivative with respect to 
X. All that remains is to calculate these expressions using eq. (22) and to 
expand all the g^s around the mid-point g„. This is a straight forward though 
tedious calculation. In this paper we will illustrate the general procedure for 
calculating S* by explicitly giving its expansion to order e^: 



go{qn ; cat) ^ goiqn,-^] + (^n 



4^ 



_ Cat 



7^1 Uri! TT + -^90 qn] TT 



+4 



16^ 



2 J 32 



32^ 
1 



+ 



,(4) 



2 J 2048^ 



5o ( ?n; Y j + 



3 /^^ . fiv\ 
2 J 



128 



^1 9^ 



^i(gn ; Cat) = (^gn; + ^c/o y) + 



8 



i/2 (^^n; 

1 



+ 



1024 



I qn, 



64 



^1 [qn, 



eAfV 
2 A 



5'2(?n ; Cat) = —5-2 ( qn, — 



(4) _ ^Af 



6144^ 



1 

128' 



^1 1 qn, ^ 



(24) 



In the above relations we expanded go up to e^, gi up to e^, etc. We also 
disregarded all the higher g^s. The reason for this is that the short time prop- 
agation of any theory satisfies oc while the gk term enters the action 
multiphed by S^'^. In general, if we wish to expand the effective action to 
we need to evaluate only gQ (up to e^^) and the remaining p — 1 functions 
gk (up to e^^"'^). The task of calculating the effective action to large powers 
of is time-consuming and is best done with the help of a standard pack- 
age for algebraic calculations such as Mathematica. Using Mathematica we 
determined the corresponding expressions for p < 9. 

Although the above system of recursive relations is non-linear, it is in fact 
quite easy to solve if we remember that the system itself was derived via 
an expansion in e^. Having this in mind we first write all the functions as 
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expansions in powers of Cat that are appropriate to the level p we are working 
at. For p = 3, we have 



5'o(gn; Cat) =1^(gn) + eNRiiqn) + e%R2{qn) 

QiiQu, ejv) = Rsiqn) + ^NRi{qn) (25) 

g2{qn]^N)^Rb{qn) ■ 

Putting this into the Eq. (24) we determine the functions Ri to R^ in terms 
of V. The p — ?> level solution equals 



go^V ^ eN— + e% 



Yl 
12 

^^^24+^^480 



92 



Yl 

24 
y(4) 

1920 



^ 240 



24 



(26) 



Note that W* depends only on the initial potential V and its derivatives (as 
well as on ejsf)- One can similarly calculate the effective action S* to any 
desired level p. We denote the p level truncation of the effective action as S^^\ 
S^^ has the property that its iV-fold amplitudes deviate from the continuum 
expressions as O(e^) 

A{a,b;T)^Af{a,b;T) + 0{e^^) . (27) 



Comparing this to eq. (6) we see that we have just derived the generalization of 
the Euler summation formula to path integrals. Just as with the ordinary Euler 
formula it gives the relation between path integrals and their discretizations 
to any given precision. 

It is important to note that one solves for the effective action at level p but once 
for all theories, i.e. the solution that is found holds for all initial potentials. 
The only requirement for the level p solution is that the starting potential 
is diffcrentiable 2p — 2 times. Solutions for larger values of p are a bit more 
cumbersome, however, they are just as easy to use in simulations. We have 
found that the growth in complexity of the effective actions with increasing p 
has little effect on computation time for p < 4:, while simulations with p — 9 
arc roughly ten times slower due to this. However, this is an extremely small 
price to pay for a gain of eight orders of magnitude in the speed of convergence. 
Expressions for effective actions up to p = 9 can be found on our web site [16]. 
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Fig. 1. Deviations from the continuum limit — A\ as functions of N for 

p = 1,2,4 and 6 for an anharmonic oscillator with quartic coupling A = 10, time of 
propagation T = 1 from a = to 6 = 1. Nmc was 9.2 • 10^ for p = 1, 2, 9.2 • 10^° 
for p = 4, and 3.68 • 10^^ for p = 6. Dashed lines correspond to appropriate 1/A^ 
polynomial fits to the data. Solid lines give the leading behavior. The level p 
curve has a l/N^ leading behavior. 



The analytical derivations presented work equally well in both the Euclidean 
and Minkowski formalism (with appropriate ie regularization) , i.e. they are 
directly applicable to quantum systems as well as to statistical ones. However, 
the Monte Carlo simulations used to numerically document our analytical 
results necessarily needed to be done in the Euclidean formalism. We analyzed 
in detail several models: the anharmonic oscillator with quartic coupling V = 
I + ^ and a particle moving in a modified Poschl- Teller potential over a 
wide range of parameters. In all cases we found agreement with eq. (27). Fig. 1 
illustrates this behavior in the case of an anharmonic oscillator. We see that 
the p level data indeed differs from the continuum amplitudes as a polynomial 
starting with l/N^. The deviations from the continuum limit |y4j^^— y4| become 
exceedingly small for larger values of p making it necessary to use ever larger 
values of Nmc so that the MC statistical error does not mask these extremely 
small deviations. For p = 6 we see that although we used an extremely large 
number of MC samples {Nmc = 3.68 ■ 10^^) the statistical errors become of 
the same order as the deviations already at ^ 8. For p = 9 this is the case 
even for N = 2, i.e. we already get the continuum limit within a MC error of 
around 10~®. 
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To conclude, we have presented an algorithm that leads to significant speedup 
of numerical procedures for calculating path integrals. The increase in speed 
results from new analytical input that comes from a systematic investigation 
of the relation between discretizations of different coarseness and that leads 
to a generalization of the Euler summation formula to path integrals. We have 
presented an explicit procedure for obtaining a set of effective actions S^^ that 
have the same continuum limit as the starting action S, but which approach 
that limit ever faster. Amplitudes calculated using the A^-point discretized 
effective action S^'^ satisfy A^p{a, b; T) = A{a, b; T) + 0{1/Np). We have ob- 
tained and analyzed the effective actions for p < 9 and have documented the 
speedup up to by conducting Monte Carlo simulations of several differ- 

ent models. Several interesting properties of this procedure follow from the 
fact that the solutions were obtained using an asymptotic expansion. These 
additional properties will be presented and discussed in a future publication. 
Extension to d > 1 is also in progress. The derivation of higher dimensional 
analogues of integral equation (21) docs not seem to present a problem. The 
asymptotic expansion used to solve it is also directly generalizable. However, 
the algebraic recursive relations that determine 5**^^^ will be more complex and 
may practically limit us to smaller values of p. 
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